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Abstract 

A  class  of  useful  difference  approximations  to  the 
full  non-linear  Navier-Stokes  equations  is  analyzed;  the 
convergence  of  the  solutions  to  the  solutions  of  the 
corresponding  differential  equations  is  established  and 
the  rate  of  convergence  is  estimated. 


ii 


Introduction.   The  Navler-Stokes  equations,  describing  the 
motion  of  a  viscous  Incompressible  fluid,  can  be  written  In  the 
dlmenslonless  form 

(1)  S^v  +  grad  p  =  -^j^j^i  +  ^^Z  +  E,  (v^  =  3^   ^j  } 

(2)  dlv  V  =  0 

where  the  vector  v,  with  components  v.,  1  =  1,2,5,  Is  the  velo- 
city, p  Is  the  pressure,  E  Is  the  external  force,  3,  denotes 
differentiation  with  respect  to  the  time  t  and  5.  denotes  dif- 
ferentiation with  respect  to  the  space  variable  x.,  1  =  1,2,3. 
Vector  quantities  are  underlined  and  the  summation  convention 
applies  to  the  Index  j. 

When  a  solution  of  these  equations  Is  required  In  some 
bounded  domain  Q,  with  boundary  ^f2,  use  Is  generally  made  of  an 
appropriate  difference  approximation.  A  new  class  of  such  ap- 
proximations was  Introduced  and  utilized  In  [1]  and  [2];  It  Is 
the  purpose  of  this  paper  to  establish  the  convergence  of  the 
solutions  of  such  approximations  to  the  solutions  of  equations 
(1)  and  (2)  In  Q. 

To  our  knowledge,  the  first  convergence  proof  for  a  dif- 
ference approximation  to  the  complete  system  (l)  and  (2)  was 
given  by  Krzhlvltskl  and  Ladyzhenskaya  (see  e.g.  [3])-   Their 
proof  gives  both  more  and  less  than  the  numerical  analyst  re- 
quires.  It  gives  more  because  it  actually  establishes  the 
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existence  of  a  certain  weak  solution  of  the  equations.   It  gives 
less  because  It  provides  no  estimate  of  the  error  and  because 
It  applies  to  a  scheme  which  Is  not  readily  applied  In  practical 
calculation.   Proofs  related  to  that  of  Krzhlvltskl  and 
Ladyzhenskaya  have  been  given  by  Temam  [4],  [5],  for  schemes 
which  are  as  yet  untested  In  practice. 

In  the  present  paper  we  shall  adopt  a  different  point  of 
view.   We  shall  assume  that  the  differential  equations  have  a 
solution  with  a  certain  number  of  continuous  derivatives. 
Armed  with  this  knowledge,  we  shall  study  difference  schemes 
which  are  rot  merely  usable,  but  even  efficient.   The  methods 
analyzed  are  based  on  the  following  observations:  Equation  (l) 
can  be  written  In  the  form 


( 1 ' )  S , V  +  grad  p  =  ^v  +  E 


where  tv  =   -v.B.v.  +  V  v  Is  a  functional  of  v;  equation  (2) 

^™  J   J   "*■ 

can  be  differentiated  to  yield 


(2')  dlv  (h^v)   =   0 

(1')  can  therefore  be  written  In  the  form 
(3)  \v  =  !P{4z  +  E) 

where  Xls  an  orthogonal  projection  operator  which  projects 
vectors  In  LpC^)  onto  the  subspace  of  vectors  with  zero 

-2- 


divergence  In  Q   and  satisfying  an  appropriate  boundary  condition 
on  ho,    (see  e.g.  [6],  [7]).   On  the  basis  of  these  remarks  the 
following  procedure  is  followed:   The  time  t  is  discretized; 
at  every  time  level  <¥-v,  then  J^(^v  +  E)  are  evaluated;  this 
yields  an  approximation  to  ^, v  which  is  used  to  obtain  v  at  the 
next  time  level. 

As  will  become  apparent  in  the  course  of  this  work,  the 
author  has  not  obtained  results  as  gereral  as  he  may  have 
wished.   A  convergence  proof  in  both  the  maximum  and  Lp  norms, 
with  a  suitable  error  estimate,  has  been  obtained  only  for  the 
special  problems  in  which  the  boundary  conditions  are  replaced 
by  periodicity  conditions.   This  proof  is  presented  in  the  next 
two  sections;  first  the  discrete  analogues  of  the  operators 
grad,  div  and  y  are   described  and  studied;  these  operators  are 
then  used  to  present  and  analyze  a  difference  scheme  for  the 
periodic  initial  value  problem.   The  mixed  initial  value -boundary 
value  problem  is  briefly  discussed  in  a  final  section. 


Preliminaries;  the  operators  D,G  and  P.   We  assume  that  equa- 


tions (1)  and  (2)  have  a  solution  v,  p,  periodic  in  all  spatial 

directions;  without  loss  of  generality  the  periods  can  be 

taken  equal  to  1.   Let  £   be  the  number  of  space  dimensions;  Q, 

is  then  the  cube  0<x.  <  1,  i  =  !,...,£.      We  cover  fi  by  a 

rectangular  grid  and  assume  that  the  mesh  widths  in  all  directions 

are  equal  to  the  same  small  number  h.   The  set  of  all  mesh  nodes 

is  denoted  by  Q,.;   0,  is  the  set  of  nodes  in  the  interior  of  fi 
"^   h   h 
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and  ^il  Is  the  set  of  nodes  on  the  bo-undary  of  fi.   n°  +  hO.     =  Q,    . 

N  =  h   +  1  Is  the  number  of  mesh  points  In  each  space  direction. 

Let  f  be  a  scalar  function  and  let  u  be  a  vector  function 

with  components  u.,  defined  at  the  points  of  fi,  .   Let  z  =  [q,r] 

(if  i  =  2)  or  2  =  [q,r,s]  (if  i  =  3)  be  a  point  in  Q     with 

coordinates  qh,  rh  (qh,rh,sh  if  ^  =  3) •   The  values  of  f,  u. 

at  z  are  denoted  by  f  ,  u. /  %  or  f    ,  u^  /    n  (f     >  u./     \, 

^      z   l(z)     q,r'   l(q,r)  ^  q,r,s'   i(q,r,s)' 

if  i  =  5)'   The  periodicity  conditioriS  become 


f        =  f 
q+N-l,r    q,r 


etc . 

The  inner  product  Is  defined  for  scalar  functions  f,g,  by 


and  for  vectors  u,  v  by 


(u.v)  =  YH    ("I'V^) 


1=1 

As  usual,  we  set 


||f|l  =/(f7fT  ,     \\v\\  =  /Xu7uJ 


The  shift  operators  S+.  are  defined  by 


-4- 


^+l^q,r  -  ^q+l,r 


^I2^q,r  -  ^q,r+l 


with  similar  definitions  in  the  three  dimensional  case.   The 
difference  operators  D^,  D_^,  D  .  are  defined  by 

D+1  =  (S^.1  -  l)/h 
D_i  =  (I  -  S_^)/h 
^01  =  ^^+1  +  ^-1^/2  =  (S^.  -  S_^)/2h 

where  I  is  the  identity.  ^,^>  D_.  and  D  .  are  respectively  the 
forward,  backward  and  centered  difference  operators  in  the  1-th 
direction. 

Let  D  and  G  denote  respectively  the  discrete  approximations 
to  the  operators  dlv  and  grad.   Both  D  and  G  employ  centered 
differences,  i.e.  for  a  vector  u_  on  Q,,    we  set 


Du  =  D  .-u. 
—    oj  j 


and  for  a  scalar  function  \)   we  set 


Gi^  =  Dol^ 


With  these  definitions,  the  following  identities  can  be  readily 
verified: 
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(k)  (Du,e)  =  0 


where  e  =  1  at  all  points  of  n,  ,  and 


(5)  (Du,(t))  +  (u,G(t))  =  0 

for  all  u  and  \).      These  are  the  analogues  of  the  identities 


/  dlv  u  dx  =  0 
Q, 


and 


/  (})  dlv  u  dx  +  /  u  grad  \>   dx 
=  /  dlv  (j^u  dx  =  0 


which  hold  for  smooth  periodic  functions  u  and  (j)  on  Q.   For  u,)}) 
periodic  and  three  times  continuously  differentiable .  we  have 


\\g\>   -   grad  ^\\    =  O(h^) 


||Du  -  dlv  u||  =  O(h^) 


We  shall  now  discuss  some  consequences  of  our  systematic  use  of 
centered  differences.   Let  t^  be  a  function  on  Q^,    let  z_.  z^ 
be  two  points  a  distance  2h  (modulo  1)  apart,  and  let  z^  be  the 
point  on  the  line  joining  z  and  z   and  at  a  distance  h  from 
each.   One  of  the  components  of  Gi   at  z^  is  a  linear  combination 
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of  lA   and  1^   ;  we  describe  this  situation  by  saying  the  z  and 

Z  2  "T 

+ 
z_  are  G-connected.   We  say  that  points  z,z'  belong  to  the  same 

G-chaln  If  there  exist  points  z, ,Zp,...,z   such  that  any  two 

successive  points  In  the  sequence 


Z^Z-|  jZp,  .  .  .  ,ZjZ 


are  G-connected.   Clearly  fi.  Is  the  union  of  some  number  L  of 
disjoint  G-chalns.   If  N  Is  even  L  =  1;  If  N  Is  odd  L  =  2"^. 
(Had  we  allowed  the  numbers  of  mesh  points  In  the  several  space 
directions  to  differ  from  each  other,  we  would  have  found  that 
L  =  2""",  1  between  0  and  H) .      The  following  facts  can  now  be 
verified:  (l)  When  Gif)  Is  given,  ^   Is  determined  only  up  to  L 
arbitrary  constants;  and  (11)  The  sum  on  the  left  hand  side  of 
the  Identity  (4)  can  be  separated  Into  L  partial  sums,  each 
vanishing  separately. 

We  are  now  ready  to  prove  the  following  discrete  analogue 
of  a  well-known  decomposition  theorem: 

Theorem  1.  Let  u  be  a  vector  on  a  satisfying  the  periodicity 
conditions;  then  there  exist  a  unique  periodic  vector  u  and  a 
periodic  function  ^   such  that 


(6)  Du-^  =  0 


D 


(7)  u  =  H  +  # 


at  all  points  of  n,  ,  with 
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(8)  {v^,G\>)    =   0 


Proof;   If  u  satisfying  equation  (6)  exists,  then  equation  (8) 
Is  clearly  satisfied  because 

(u^,G^)  =  -(d/,4))  =  0 

We  already  know  that  ^   in  equation  (7)  can  be  determined 
only  up  to  L  arbitrary  constants.   To  lift  this  indeterminacy 
we  can  Impose  L  additional  conditions;  for  example,  we  can  number 
the  G-chalns  and  require  thar 

(9;  5   3_«l>^  =  0     1  =  1,...,L 

where  )   .  denotes  summation  over  the  i-th  G-chain. 

The  theorem  is  proved  by  verification  of  the  Predholm 
alternative.   Let  q  be  the  number  of  points  in  Q.,    and  q^  the 
number  of  points  on  bO,^    (q^  +  q^  =  N  ) .   There  are  q^  +  qy  2 
values  of  ^   and  i(q  +  q^2)  components  of  u  to  determine. 
Equation  (6)  represents  (q  +  qy'2)  equations  related  by  the  L 
identities  (4),  i.e.  q^  +  qV^  -  L  independent  equations. 
Equation  (?)  represents  ^(q  +  <1^2)  relations;  together  with 
equations  (9)  the  number  of  equations  equals  the  number  of 
unknowns. 

Squaring  (7)  and  using  (8)  we  obtain 


(10)  ||u||2.  ||u^ll2+  ||G^||2 
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Therefore,  If  u  =  0,  then  u  =  0,  Gtj)  =  0  and  (()  =  0;  this  proves 

the  theorem. 

Let  H  be  the  space  of  periodic  vectors  defined  on  a  ,  let 

H|^  be  the  subspace  of  periodic  vectors  v  such  that  Dv  =  0,  and 

let  H„  be  the  subspace  of  vectors  of  the  form  Gh,    h   periodic; 
u  — 

Theorem  1  states  that  H^  and  H,^  are  orthogonal  to  each  other, 
and  that  H  Is  their  direct  sum.   Let  P  be  the  orthogonal  projeC' 
tlon  projecting  H  on  H^^;  (?)  can  be  written  In  the  form 


(11)  u  =  Pu  +  G(t) 
We  obviously  have  for  all  u 

(12)  l|Bill<l|u|| 

P  Is  the  discrete  analogue  of  j.      (see  equation  (5))-   Given 
u,  It  Is  a  fairly  simple  matter  to  evaluate  Pu;  efficient 
methods  for  so  doing  were  described  In  [2].   For  the  sake  of 
completeness,  we  exhibit  a  method  for  finding  Pu  which,  albeit 
inefficient,  has  the  merit  of  simplicity.   Consider  the  Iteration 
system 

w"^l  =  u  -  G^       m  >  1 
^^^   =  ^"^   -   eDw"^l   m  >  1 

where  w"^,  \>^ ,   m  >  1,  are  periodic  and  9  is  a  parameter.   It  is 
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2/«2 


readily  verified  that  for  0  <  h  /i  ,  w  converges  to  Pu_  as  m 
increases  for  all  initial  guesses  w  ,(|)  . 

To  conclude  this  section,  we  prove  a  number  of  inequalities 
which  will  be  needed  ir  later  sections.    We  start  with  a  dis- 
crete analogue  of  the  Poincare  inequality.   Consider  the  case 
&  =  2.      Let  1^  be  a  function  defined  on  Q,,,    and  let  z  =  [p,q], 
z'  =  [pSq']  be  two  points  on  the  i-th  G-chain;  p'  =  p  +  2m,, 
q '  =  q  +  2m2 .   We  have , 


m-j  -1 


P'<^^&  (^l^Vl+2k,q  •  2h 
m2-l 
+  g^  (G2^)p.,q+l+2k  • 


2h 


Therefore 


■f   ,      ,    -  if       I  <  4 
P',q'    P,q   - 


<  8 


N-1 


N-1 


g:|Gl^lk,qh  +  |=l^2^lp-,lc 


N-1 


N-1 


ferr   1    ^'^      fa:   ^    p  '^ 


h 


h 


where  the  relation  (N-l)h  =  1  Is  used.   We  multiply  both  sides 
by  h  and  sum  over  all  [p,q]  and  [p',q']  in  the  i-th  G-chain, 


giving  points  on  dfi,  the  weight  ^,    and  obtal 


n 


2   ,  2 


Tzyj:z,\%,^^'  ^  zZih^iZii^i;,,q.h^  +  ^yz,^^,^^^%^. 


q 


,h' 


<3 


<  87~^h^||GT/^||  <  8||G-^ir 
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where  >    denotes  summation  over  the  1-th  G-chaln. 

Let  N^  be  the  number  of  points  In  the  1-th  G-chaln.   We  have 


2  1.  1  .t2  1 


N,  >  {N-l)'^>^N-^ 


Therefore 


^1^'  >  ^ 


^ 


and 


E^Ill^ip.qh^  <  2(IZi*p,/)^  +  8||0*|P 


Summing  over  all  G-chalns  we  obtain 


.2x2 


(12)  |Wr<2L     YZ        (ZZA   a^    )      +   ClII^^H 

G-chalns  ^  ^'^  ^    ~ 

o 
where  G-,  =  8l  .   A  similar  Inequality  can  be  derived  In  the 

three  dimensional  case,  with  C-,  =  12L  .   The  Inequality  (12) 

can  now  be  used  to  prove  the  following  theorem: 

Theorem  2.   Let  u  be  a  periodic  vector  defined  on  Ct,.      Then  the 
—  h 

following  Inequality  holds 


(15)  IIu-phII  <  C-lI|Du_ 


where  C-,  Is  a  constant  Independent  of  u_  and  h 
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Proof:   By  Theorem  1,  u  -  Pu  =  G(t)  Is  In  H^.   Let  Gi'   be  an  arbitrary 
unit  vector  In  H^,  (IIGt^H  =  l).   V'  Is  determined  only  up  to  L 
arbitrary  constants  which  can  be  chosen  so  that 

IZi^z  =  0     1  =  1,...,L 
We  have 

{G\>,Gi')    =    (u    -   Pu,Gt)    =    (u,G'^)    =    -(Du,t) 
Hence,    using    (12),   we   obtain 

|(G^,G^)|    =    |(Du,^)|    <    llDul!    |1^||<   C^llDul! 


Since  Gif   Is  an  arbitrary  unit  vector  in  H_ ,  (13)  follows 

G 


Solution  of  the  periodic  Inltlal-value  problem.   In  this  section 
a  scheme  for  finding  periodic  solutions  of  equations  (l)  and 
(2)  will  be  analyzed.  The  particular  scheme  discussed  has  been 
singled  out  because  It  resembles  schemes  the  author  has  used  In 
actual  computation  (see  [2]);  It  will  be  evident  that  the  analysis 
applies  to  wide  classes  of  schemes.   We  shall  again  simplify 
notations  by  writing  the  equations  for  the  two  dimensional  case; 
the  scheme  as  well  as  the  proofs  generalize  to  the  three  dimen- 
sional case  without  further  ado. 

Let  u,  with  components  u.  ,  be  the  computed  velocity,  let  ir 
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be  the  computed  pressure,  and  let  k  be  the  time  step.   We  write 
u^  =  u(nk),    IT     =  Tr(nk) ,    etc. 


At  the  time  t  =  0  a  periodic  velocity  field  u  Is  assumed  given. 
(More  will  be  said  later  about  the  proper  choice  of  u°).   Given 
u",  u^^  Is  evaluated  In  three  fractional  steps: 


(14a)    uf  1/^  =  uj  -  K^ol^f  ^^^  +  ^+l^-l"f  ^""^ 
(14b)   uf  2/^  =  uf  2/3  .  ku^D^^^f  2/3  ^  kD^^D_2uf  2/3 
(14c)     u"+l  =  P(u^+2/3  ^  i^E^+l) 

with  u"^^/^,  u"'^^/^  periodic. 

Equations  (l4a)  and  (l4b)  can  be  rewritten  In  the  form 

(15a)  (I  -  kQ3^(u"))u"+V3.^n 

(15b)  (I  -  kQ2(u"))u^^^/^=u"+l/^ 

where  Q  (u^),  Qp(2i")»  ^^®  linear  operators  dependent  on  the 
parameters  u^/  \.   Equation  (l4c)  can  be  rewritten  In  the  form 

/-.r-    \  n+1  ,  ,  _,  n+1    n+2/5  ,  i -r-n+l    /y.   n+1    ^>, 

(15c)     u    +  kGir    =  u   '   +  kE      (Du    =■  0) 

which  defines  ir^    ,    the  computed  pressure  at  the  time  (n+l)k. 
(yn+2/:5  corresponds  to  u^""^  In  the  notations  of  [1]  and  [2]). 
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It  can  be  seen  that  the  vector  (u^       H  )/^  approximates  ^u 

and  that  equation  (l4c).  which  is  equivalent  to 

(u^+l  -  u^)/k  =  p{(h"-'^/^  -  u")/k  +  e} 

Is  the  discrete  analogue  of  equation  (3). 

The  task  now  at  hand  Is  to  prove  that  u^"^-"-/^,  u""*"^^^,  u"^-*" 
exist.  I.e.  that  the  operators  (I  -  kQ. (u  ))  are  Invertlble  when 
u  Is  chosen  appropriately,  and  that  the  vectors  u^  converge  to 
the  solution  v(nk)  of  equations  (l)  and  (2).   We  start  by  showing 
that  equations  (l4)  are  consistent  with  equations  (1)  and  (2); 
this  Is  the  content  of  the  following  lemma: 

Lemma  1 .   Let  k  =  0(h  ),  and  assume  that  equations  (1)  and  (2) 
have  a  periodic  solution  v,p,  which  has  continuous  derivatives 
up  to  order  five   in  the  interval    0  <  t  <  T.   Then  there 
exist  two  times  continuously  dlfferentlable  vectors  w  ,  w   '    , 
^n+2/3  (0  <  nk  <  T)  such  that 

(l6a)         (I  -  kQ^(w^))w"+^/^  =  w"  +  O(k^) 

(16b)       (I  -  kQ2(w"))w""'^/^  =  w^-^^^^  +  0(k2) 
(16c)         w"+l  =  P(w"+2/3  +  ksH+l^  ^  Q^j^2) 


with 


(17)  ||v"-  w^||=  0(k) 
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Proof:   We  simply  construct  the  required  functions.   We  have 


0        2 
Dv  =  dlv  V  +  >   ^-j-  OqV„  +  0(h  ) 


Therefore,  puttlrf 


h^   2  n 
(l8)  Wg  =  Vg  -  xr  '^ft'^ft   P  =  1,...,Z      (no  summation  over  p) 


we  obtain 

Dw""  =  OCh"^)  =  0(k2) 
Equation  (17)  Is  clearly  satisfied,  and  by  Theorem  2 
(19)  ||w"  -  Pw"||=  0(k2) 

We  now  put 


n+1 


wfV5  =  ^n+1  _  ka^vf  l+kv^a^vf  1  +  ka.p"+l  -  kE? 

n4-2/5  ^  n+1  ^  j^.   n+1  _  ^n+1 
1        1       1         1 


Equations  (l6a)  and  (l6b)  are  clearly  satisfied,  and  since 


Gp  =  grad  p  +  0(h  )  =  grad  p  +  0(k) 


we  have 
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P(k  grad  p)  =  O(k^) 


On  the  other  hand,  It  can  always  be  assumed  that  dlv  E  =  0  (see 
e-g-  [7]) J  ard  therefore  by  Theorem  2 


I  I  kE   -  kPEll  <  C^k||DE||  =   O(k^) 


These  equations,  together  with  equation  (19)  show  that  (l6c)  Is 
satisfied,  and  the  lemma  Is  proved. 

We  shall  use  w"^  as  a  comparison  vector.  I.e.  we  shall  prove 
that  I  lu'^  -  w^^ll  Is  small,  and  use  (17)  at  the  end  of  the  argument 
to  show  that  ||u  -  v  ||  Is  small.   The  lemma  assumes  that  v  has 
continuous  derivatives  up  to  order  five.   Had  we  assumed  the 
exlstei.ce  of  only  four  continuous  derivatives,  the  error  term 
Ir  (l6c)  would  have  been  of  order  kh.   This  is  sufficient  for 
convergence;  however,  the  proof  becomes  somewhat  more  complicated 
and  we  shall  be  content  with  the  assumption  of  the  lemma.   It 
should  be  pointed  out  that  even  for  the  simplest  linear  parabolic 
equations  it  is  t ecessary  to  assume  the  existence  of  continuous 
derivatives  up  to  order  four  in  order  to  obtain  the  correct 
order  of  magnitude  of  the  truncation  error. 

We  now  ir.troduce  a  second  norm,  the  discrete  maximum  norm, 
defined  for  scalar  functions  \)   by 

and  for  vectors  u  by 
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llHllmax=  "^^^  1 1^1 


.     J.  max 


We  have 
Lemma  2 . 


h-V2 


The  proof  is  obvious  from  the  definitions.   This  lemma  is  crucial 
to  the  sequel  since,  as  we  shall  see,  it  implies  that  if  u_  con- 
verges to  V  with  sufficient  accuracy  ir  the  Lp  norm,  then  u_ 
also  converges  to  v  in  the  maximum  norm. 

Lemma  3.   Let  ||u^||  ^^  <  K,  and  let  k  be  small  enough  for  the 
^^— — ^^— ^—         max  — 

inequality 


kK^/4  <  1 


to  hold.   Equatior s  (l4a)  and  (l^b)  can  then  be  uniquely  solved 

„    n+l/3   n+2/3 
for  u     J  H     • 

Proof;   Multiplying  (l4a)  by  u?"*"  '-^   we  obtain 
||.5+VJ||^  ,(„n+l/5_„nj,^^,n+l/5) 

however,  we  have 
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,    n+l/5      T^     T^        i^+l/5\  iit^         I  i2 

(!.'?+ 1/3^    ^n      ^n+l/5)|         ||^n||  ||u?^l/3||   ||d      n+l/3| 

^1  loll  /i^iiiii  "max  "1  "   "   d   1  ' 


On    the   other  hand. 


^01  =^^^1  +  I^-l) 


therefore 


i%i^^'^'ii<  i(iii'+i"f'^^ii  +  iii'-i^r'^'ii)=  iKi^'^'ii 


and 


and  hence 


(20)  iK+l/^lld    -^)    <    \Wl\\ 


The  existence  and  uniqueness  of  u?  ^      follow  by  the  Fredholm 
alternative.   The  existence  and  uniqueness  of  u.  '      are  estab- 
lished In  the  same  way. 

Lemma  4.   Let  ||u"||  „^  <  K,  and  let  k  be  small  enough  for  the 
Inequality 


-18- 


kK^/2   <   1 


to  hold;    then,    if  k  =   0(h2)^   ^q   have 


(21)  llu"*^    -  w"+^||<    (1  +   kC2(K))|Iu''    -  w"||+   C^kh^ 


where  C^  Is  a  constant,  and  C^Ck)  is  a  constant  whose  magnitude 
depends  on  K. 

Proof:   Subtracting  (l4a)  from  (15a)  we  obtain 

n+1/3  _  n+1/3  ^  _^n        (^n+l/3  _  ^n+l/3) 
1         1  1  ol^  1         1     ^ 

,  /  n    nv^    n+l/3 

M-ultlpllcation  by  u?  ^   -  w?  ^   and  manipulations  similar  to 
those  in  the  proof  of  Lemma  3  yield 


where 


M-,  =  max  max  jS  w^"^^^] 
^   0<t<T  Q  ^ 


Similarly,  we  obtain 


■19- 


where 

Mp  =  max  max  |  S^w'^'  '     \ 
0<t<T  fi 

and  hence 

where  CoCk)  depends  on  K.   Finally 

u"+l  -  w"+l  =  P(u"+2/3  _  ^n+2/3)  ^  o(i,2) 

and,  therefore,  using  (12) 


l^n+1    _  ^n+l||^    ii^n+2/3    _  ^n+2/3||^   o^i^2^ 


<    (1  +   kC2(K))||u"   -  w^^ll  +  C^kh^ 


and  the  lemma  is  proved. 

Let  u°  be  the  initial  value  of  u,  for  use  in  equations  (l4) 
We  assume  that 


(22)  ||u°-  w°||=  C^h2 


This  can  be  achieved  for  example  by  putting 
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o    o 

U   =  V 


Let  W  be  defined  by 


W  =  max  max  max  |w. 
1   0<t<T  fi    ^ 


Let  C  be  the  largest  of 


C^,  Ci^  and  C^{2]f]) 


Assume  k  Is  so  small  that 


kC  «  1 


and  that  h  is  smaller  than  h  ,  where 

o 


{c,i)hl^-^^^^  =  ifeWj      e<  1 


max    ,_, 


Equation    (22)   and  Lemma  2   then   show  that 


II"°llmax  <  W+^W<  2W      . 


By  Lemma  3  u     exists  and  by  Lemma  4  we  have 


||u     -  w    II  <    (1  +  Ck)Ch     +  Ckh' 


Therefore 
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1 1   1 

u  -  w 


lax  ^  ¥^   +  ^^^^  +  f  kW  <  ^W  <  W,    Ih'^IUx  <  ^^  ' 


2 
and  we  can  evaluate  u  .   In  general  we  have 


(23)    ||u''"^^-w"+^||<    (l+Ck)"+^Ch^  +    [1  +    (l+Ck)+...+  (l+Ck)"]Ckh^ 


<  2eC("+l)^  max    (C,l)h2 


and 


(24)  ||u"+l   -  w"+l||^g_^  <     6We^^    (t  =    (n+l)k) 


Let   T     be   defined  by 


e-o  .  1 


and  let  T^  =  min  (T,T^).   Inequality  (24)  shows  that  for 


0  <  t  <  T^ 


IML«.  <  2W 


max 


and  hence  for  0  <  nk  <  T-,  ,  u""^   exists  and 


(25a)  ||vP+^    -  w'^+^lll  2  max    (G,l)e^V 


as  well  as 
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(25b)     l|u«l  -  w-^lll  „_  <  2   max  (cDe^hC-^'/a 


max  — 


If  T-,  <  T,  I.e.  If  the  Inequalities  (25)  hold  for  a  time 
Interval  shorter  than  the  time  Interval  for  which  the  solution 
of  the  differential  equations  has  5  bounded  derivatives  and  for 
which  a  numerical  solution  Is  required,  the  above  process  can 
be  restarted  at  t  =  T^ ,  to  yield  convergence  for  the  whole 
finite  Interval  0  <  t  <  T. 

Bearing  In  mind  the  definition  of  w  and  equation  (17) ,  we 
obtain  the  following  theorem: 

Theorem  3-   Let  equations  (1)  and  (2)  have  a  periodic  solution 
with  continuous  derivatives  up  to  order  five  for  0  <  t  <  T.   Let 

k=0(h  );  if  |[u°-w°||,  k  and  h  are  sufficiently  small,  equations  (l4) 
have  a  unique  solution  which  converges  to  the  solution  of  (l) 

and  (2)  In  both  the  Lp  and  maximum  norms.   The  error  in  the  Lp 

2 
norm  Is  of  order  h  ;  the  error  In  the  maximum  norm  is  bounded 

by  0(h)  In  the  two  dimensional  case  and  by  0(/h)  In  the  three 

dimensional  case. 

Theorem  3  and  Its  proof  can  be  summarized  as  follows:   Let 

u",oo"  be  vector  functions  defined  for  z  In  fi,  and  for  n  such 
— z  — z  h 

that  0  <  nk  <  T-,  ;  Introduce  the  "space -time"  maximum  and  'L,^   norms 


1^1  max  T  =   "^^^   "^"11  max 
max,i^   0<nk<T-j^     ^^^ 


lyL     _       max      ||u"| 
•^1   ~   Q<nk<T^ 


-25- 


The  equations 


u   given  , 


define  a  mapping  oj  ->  u.   This  mapping  maps  the  maximum  norm 
sphere 


l^lmax,T^  <  2||vlIUx,T^ 


Into  the  Lp  norm  sphere 


V2 
1~   -  ..^..,.-^' 


IIh-w|It  <  IlwlUx,T  h 


For  I  |u  -  w  ||,k  and  h  sufficiently  small,  this  mapping  has  a 
unique  fixed  point  which  Is  a  solution  of  (l4)  and  lies  close 
to  V,  the  solution  of  (l)  and  (2). 

In  our  analysis  we  have  neglected  the  effect  of  round-off 
error  and  of  the  errors  arising  from  the  possibly  Incomplete 
Iterative  evaluation  of  u'^"^  =  P(u'^"''  '     +  V^^      ).      It  Is  obvious 
however  that  the  aiialysls  remains  valid  if  the  round-off  errors 
are  of  order  k  and  provided  u    is  approximated  by  a  vector 
(u    )   such  that 
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D(h"^1)*  =  0(k2) 


Furthermore,  In  the  dlmensionless  variables  used  In  this  paper  the 
effect  of  the  Reynolds  number  R  on  the  error  is  not  In  evidence. 
Clearly  C  depends  on  R  and  Increases  as  R  increases;  I.e.  as  R 
increases  k  and  h  have  to  be  reduced  for  accuracy  to  be  preserved. 
Finally,  It  Is  clear  that  the  results  of  this  section  apply  to 
certain  other  quasi-linear  equations  besides  the  Navler-Stokes 
equations,  provided  the  boundary  conditions  are  homogeneous. 
In  this  sense,  our  results  generalize  the  work  of  M.  Lees  (see 
e.g.  [  8]),  who  considered  equations  with  nonlinear  terms  of  a 
simpler  nature. 


The  mixed  initial  value -boundary  value  problem.   The  main  Inter- 
est of  methods  such  as  those  considered  in  this  paper  lies  in 
their  applicability  to  mixed  Initial  value -boundary  value  prob- 
lems.  Schemes  similar  to  (l4)  have  been  successfully  applied 
by  the  author  to  a  variety  of  such  problems  (see.  e.g.  [2]). 
The  convergence  proof  however,  becomes  considerably  more  diffi- 
cult in  the  presence  of  boundaries. 

Consider  in  particular  the  problem  of  solving  equations  (l) 
and  (2)  in  a  domain  Q.,   with  v  given  and  with  the  boundary 
condition 

(26)  V  =  0   on  Sfi  . 
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Operators  D  and  G  can  be  constructed  so  that  the  Identities  (  ^  ) 
and  (  5  )  are  satisfied  and  Theorems  1  and  2  hold.   D  and  G  thus 
constructed  employ  centered  differences  except  on  bQ.      On  hO, 
one-sided  first  order  differences  are  used  whenever  the  use  of 
centered  differences  would  require  functional  values  at  points 
outside  n.   The  projection  P  associated  with  G_   and  D  Is  ortho- 
gonal In  the  space  of  functions  satisfying  (26).   The  proofs  of 
all  these  statements  take  Into  account  the  fact  that  the  number 
of  G-chalns  is  2   independently  of  the  number  of  points  in  the 

mesh. 

Difficulties  arise,  however,  when  one  approaches  the  con- 
vergence proof  proper.   It  is  clear  from  the  proof  of  Lemma  1 
of  the  last  section  that,  were  one  to  use  schemes  such  as  (l4)  , 
one  would  have  to  impose  on  u^"^  /^,  u"^  ^^  inhomogeneous 
boundary  conditions  of  the  form 


(27a)  u^-*-^^  =  kG-rr^   on  hQ^ 


(27b) 


u^+2/5  ^  i^Q^n   ^^  ^^^ 


where  tt  is  the  pressure  computed  at  time  t  =  nk.   Such  a  pro- 
cedure has  indeed  been  followed  in  practice.  Unfortunately,  in 
the  presence  of  inhomogeneous  boundary  conditions  the  author  has 
not  been  able  to  establish  the  analogues  of  Lemmas  3  and  4. 
Moreover,  the  construction  of  w  in  Lemma  1  does  not  carry  over 
to  the  present  problem,  since  w,  as  given  in  the  last  section, 
does  not  satisfy  the  imposed  boundary  conditions.   Both  difficulties 
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stem  from  the  fact  that  In  the  presence  of  boundaries  the 
operators  5*  and  V   (llnpar  part  of  2  )  do  not  commute,  nor  do 
their  discrete  analogues.   This  Is  reminiscent  of  other  situa- 
tions In  numerical  analysis  where  the  non-commutatlvlty  of 
certain  operators  hinders  the  analysis  of  fractional-step 
methods  without  detracting  from  their  practical  usefulness. 

It  Is  nevertheless  possible  to  develop  schemes  for  which 
convergence  in  the  Lp  norm  can  be  established.   As  an  example, 
consider  the  following  scheme  with  two  fractional  steps: 

(28a)    uJ+V2  ,  ,n  .  ^5^j-|(s^p  ^  S.p)(.gD„pU^,^V2) 

(28b)    h""^/^  =0  on  in^ 

(280)    u"+l  =  p(«"-^/2  ^  ^^m-1) 

It  is  clear  the  homogeneous  boundary  condition  (28b)  contains 
an  error  of  order  k.   However,  since  the  number  of  mesh  points 
on  the  boundary  is  0(h)  times  the  number  of  mesh  points  in  the 
whole  domain,  some  accuracy  in  the  Lp  norm  will  be  preserved. 
We  shall  indicate  how  one  can  establish  that  in  the  Lo  norm  the 
solution  of  (28)  converges  to  the  solution  of  the  Navier-Stokes 
equations  which  satisfy  the  correct  boundary  coridltions.   u  ,  as 
given  by  (28),  therefore  assumes  the  Imposed  boundary  conditions 
in  a  weak  sense.   It  Is  clear  that  the  estimates  we  shall  derive 
will  not  do  Justice  to  the  accuracy  of  the  method. 
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One  can  verify  the  following  Identity 


(f^Ii:|(S+p+  S_p)(u^D^pf))  =  0 


which  holds  for  all  f  provided 


Du  =  0   in  Q, 

—  h 


and  u  =  0  on  the  boundary.  This  of  course  is  a  discrete  analogue 
of  the  identity 


/  fu.a.fdx  =  0 


which  holds  whenever  div  u  =  0  in  f^  and  u  =  0  on  the  boundary. 
Using  this  identity  the  following  inequalities  can  be  established; 


u 


"+1/2,1  <  „^n, 


and 


n+1 
1=0 

If  we  assume  that  equations  (l)  and  (2)  have  a  solution  v 
with  continuous  derivatives  up  to  order  four,  this  inequality 
can  be  used  to  show  that  if  k  =  0(h  ) 


(29)        ||u^  -  v"||<  constant  /h,    0  <  nk  <  T 
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For  two  dimensional  problems  one  can  replace  (28a)  by  an  explicit 
scheme  (which  does  not  require  intermediate  boundary  data 
such  as  (28b)).   For  small  enough  Reynolds  number  and  provided 
k  <  h  /4  one  can  then  derive  an  estimate  similar  to  (29). 
Furthermore,  the  scheme  (28)  can  be  modified  so  that  a  conver- 
gence proof  of  the  Krzhivltskl-Ladyzhenskaya  type  becomes  possible 

Since  neither  the  scheme  (28)  nor  its  modifications  are  of 
any  particular  practical  significance,  details  and  proofs  are 
omitted,  (it  should  be  pointed  out  however,  that  the  system  of 
linear  equations  (28a)  can  be  solved  by  successive  relaxation, 
provided  the  relaxation  factor  o)  is  sufficiently  small.   For 
proof,  see  [  9  ] ) • 

In  ending,  the  author  would  like  to  make  some  comments  on 
the  preceding  proofs.   First  of  all,  he  would  like  to  state  his 
belief  that  the  value  of  a  scheme  such  as  (l4)  lies  in  its 
practical  usefulness,  not  in  the  possible  existence  of  a  con- 
vergence proof.   The  value  of  the  convergence  proofs  lies  in  the 
fact  that  they  contribute  to  the  understanding  of  the  numerical 
processes  performed  on  the  computer. 

The  proof  of  this  paper  requires  the  existence  of  four  or 
five  continuous  derivatives  of  v  and  p.   Furthermore,  the  error 
Increases  as  the  bounds  on  the  required  derivatives  Increase. 
This  situation  is  inherent  in  the  very  nature  of  difference 
schemes;  as  a  result,  it  is  highly  Improbable  that  a  flow  con- 
taining a  strong  cascade  process,  i.e.  a  process  in  which  energy 
is  transferred  from  large  to  small  eddies,  can  be  adequately 
described  by  a  difference  method,  for  indeed,  such  flows  are 
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characterized  by  a  rapid  increase  in  the  higher  derivatives. 
This  of  covrse  excludes  turbulence  from  the  range  of  application 
of  difference  methods. 

Pirally,  it  has  been  claimed  by  several  authors  that  the 
non-linear  terms  in  the  Navier-Stokes  equations  must  be  cast 
in  "conservation  law"  form,  i.e.  in  a  form  which,  were  it  not 
for  the  presence  of  viscous  dissipation,  would  exactly  conserve 
certain  global  properties  of  the  flow.   The  author  knows  of  no 
good  reason  for  following  this  procedure  and  has  not  endeavored 
to  do  so.   Conservation  laws  have  indeed  well  known  and  important 
applications;  in  the  present  problem  those  special  forms  of  the 
non-linear  terms  may  mask  a  breakdown  in  the  validity  of  the 
schemes  and  create  an  illusion  of  accuracy  without  the  reality. 
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